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ABSTRACT 



Computer simulations were performed to locate the equilibrium 
positions and binding energies of interstitial He, Ne, Ar , Kr , and 
Xe atoms in a tungsten crystal. Heavy interstitial atoms in tung- 
sten share a lattice site with the atom that normally occupies that 
site and form what is called a split interstitial. Three charac- 
teristic interstitial sites were located relative to each lattice 
site tested. The distance of the impurity atom from the site was 
seen to vary roughly inversely with its mass, and the displacement 
of the lattice atom increased with the mass of the impurity atom. 

The foreign atom in its interstitial position was tested to 
determine the minimum initial kinetic energy needed to escape the 
lattice, as well as the optimum escape direction. The minimum 
energy may be interpreted to be the binding energy of the defect. 

A comparison of experimental binding energies from Kornelsen and 
Sinha and simulated binding energies indicates the model gives 
realistic results . 
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I . INTRODUCTION 



The advent of modern, high speed digital computers has led to 
the application of computer simulation techniques to many differ- 
ent types of physical systems. One such application is the modeling 
of the situation that occurs when a foreign atom interacts with a 
metallic crystal lattice. In general, such modeling can be broken 
down into two basic areas, dynamic simulation and static simulation. 
As an example of the former, radiation damage has been studied by 
the simulated firing of an atom or ion onto a crystal face. Other 
examples are sputtering simulations Cl, 2 , 3] in which the incoming 
particle causes surface atoms to be ejected; and channeling simu- 
lations [4] in which the ranges of ions travelling in crystal 
lattices are calculated. Static simulations, on the other hand, 
have been concerned with the equilibrium positions in the lattice 
after point defects such as replacement atoms, interstitial atoms, 
and vacancies have been introduced. Examples of this type of 
simulation can be found in [5,6,7]. This present research uti- 
lized aspects of both static and dynamic simulation techniques. 

The goal of this research was to correlate the results of experi- 
mentally determined binding energies of point defects in a tungsten 
lattice [8,9] with the results obtained by computer simulation. 
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II. NATURE OF THE PROBLEM 



A. HISTORICAL BACKGROUND 

An investigation of radiation damage events by computer simu- 
lation techniques of crystalline behavior was published in i960 by 
Gibson, Goland, Milgram, and Vineyard, hereafter referred to as 
(GGMV) [5]- This Brookhaven National Laboratory investigation set 
forth the criteria that must be satisfied before the simulation 
method could be applied to crystals. Such factors as potential 
energy functions, forces, crystallite sizes, computation methods, 
choice of time intervals, and computer limitations were discussed. 
The crystal lattice modeled in their research was metallic copper, 
a face-center cubic (fee) structure. The potential functions used 
in the calculations was a Born-Mayer repulsive potential, with the 
necessary cohesive forces applied on the boundries of the crystal- 
lite. In integrating the equations of motion, the Brookhaven 
group used the central difference procedure. The optimum choice 
for timestep duration, At, was shown to be of critical importance 
in the integration scheme. The energy of the strongest interac- 
tion governed their choice of the above parameter. The static 
results obtained by GGMV confirmed the existance of the (lOO) 
split interstitial in the fee lattice. Their dynamic results 
described collision chains and focusing phenomena in crystallites 
struck with energetic knock-on atoms. 

Additional crystal simulation studies were performed by 

R.A. Johnson C6,10,ll]. In Ref. 6, he investigated point defects 
in a copper lattice using Born-Mayer repulsive potentials. The 
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^100) split interstitial was found to be the only stable inter- 
stitial position. He found it necessary to allow the interstitial 
to interact only with its six nearest neighbors. Atoms near the 
defect were treated as independent, while the remainder of the 
metal was treated as an elastic continuum with atoms imbedded in 
it. 

Research on body-centered cubic (bcc) crystals was undertaken 
by Erginsoy, Vineyard, and Englert (EVE) [ 7 ]. They used a com- 
posite potential function for most of their detailed work. It 
consisted of an exponentially screened Coulomb potential at small 
separations, a Born-Mayer function in the region between small 
and intermediate separations, and a Morse potential at larger 
separation distances. A split interstitial was reported for simu- 
lated bcc crystals in the (lio) direction. In their dynamic re- 
sults they reported the existence of a threshold energy for 
displacement that was highly independent of the direction of 
knock-on. Also, collisional chains in the (ill) and ( 100) direc- 
tions were found to exist. R.A. Johnson also published results 
for bcc simulations involving (1-iron and tungsten til]. The 
existence of the split interstitial in the (lio) direction was 
confirmed and crowdian migration data were discussed. The present 
investigation confirmed the (lio) split interstitial positions 
for argon, neon, krypton, and xenon. 

D.E. Harrison and associates have published several articles 
in which computer simulation of crystalline behavior has been 
investigated Cl, 2 , 12]. In a study of a fee model of copper, col- 
lision events between a copper atom and a copper lattice were 
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simulated as a function of the potential function, the enegy of 
the collision, and the location of the impact point. The inte- 
gration scheme used was an average force procedure, instead of 
the central difference procedure used by GGMV. A complete dis- 
cussion of the method can be found in Ref. 13 . A continuation of 
copper simulations was published by Harrison, Leeds, and Gay in 
1965 [ 12 ]. Another paper by Harrison and Greiling dealt with the 
ranges of heavy ions in tungsten crystals whose atoms had under- 
gone thermal displacement [ 4 l. It was found that room temperature 
thermal displacements had a negligible effect upon the collisions 
for ion energies greater than a few thousand electron volts. 
Finally, Harrison, Levy, Johnson, and Effron published results on 
computer simulation of sputtering [ 2 ]. 

Research undertaken by Vine Cl4J provided the foundation for 
this author’s present investigation. Vine used the Gay-Harrison 
model for crystal simulation as modified by Levy Cl5l, Johnson [ 16 ] 
Effron [17], and Moore Cl 8 ] . His overall objective was to cor- 
relate experimental and simulated binding energies of neon and 
argon point defects. Repulsive potential functions for neon- 
tungsten (Ne-W) and ar gon-tungs ten (Ar -W) interactions were 
used [l9]. Morse functions were not used for those interactions 
because experimental data giving the Morse parameters was based on 
homogeneous media such as tungsten-tungsten (W-W) interactions C20] 
Vine subsequently attempted to correlate results of equilibrium 
position studies for tungsten defects in a tungsten lattice using 
two different potential functions for the interstitial-lattice 
interactions. In one case, the tungsten interstitial was allowed 
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to interact using a Born-Mayer repulsive potential, and the other 
case, the tungsten defect was given a composite potential that was 
identical to that given to all the other lattice atoms. By so 
doing, he attempted to establish an empirical relationship between 
the two methods to apply to the results of neon interstitial studies 
which used only a repulsive potential. 

The results of the tungsten-tungsten interactions failed to 
provide the information needed to formulate a correction factor 
to be used in the neon-tungsten studies. In addition, the concept 
of relating the potential energy at equilibrium to the experi- 
mentally observed binding energies of Kornelsen and Sinha [9] does 
not appear to be feasible. 

B. CHOICE OF POTENTIAL FUNCTION 

The studies that were reviewed in the previous section utilized 
many different approximations to the true potential function be- 
tween atoms in a metal lattice. The problem of solving the raany- 
body interaction of a real system is approximated in the computer 
by many two-body interactions. Thus, central pairwise potential 
functions are most often used in computer simulations. GGMV 
employed a repulsive potential of the Born-Mayer form: 

V ij = ex P( A+Br ij) 

which described the repulsion of atoms at close approach. Three 
Born-Mayer potentials were investigated by GGMV to determine which 
would give the best results in their calculations. Their choice 
was one which has since been labeled the* Gibson Number Two Potential. 
In Ref. 6 , Johnson and Brown used a similar Born-Mayer potential 
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with slightly different parameters. As previously mentioned, 
another potential function that has been used in crystal modeling 
studies is the familiar Morse potential of the form: 

*ij = t>[exp [- 2 a(t i . - -2 exp { - a (t i} - r 9 )}] 

where r^ is the equilibrium distance of approach of two atoms, 
and O' and D are constants. 

Girifalco and Weizer calculated Morse parameters that would 
be appropriate for several crystal lattices [20]. In calculating 
the parameters, they attempted to express the various crystal 
properties such as cohesive energy, lattice constant, compres- 
sibility and equation of state in terms of the Morse function. 

The Morse potential constants published by Girifalco and Weizer 
have been used extensively in simulating the potential functions 
and forces in lattices of homonuclear atom systems. 

The Born-Mayer potential and the Morse potential are useful 
over specific internuclear separation distances. The Born-Mayer 
potential is useful at strongly repulsive separations, i.e. short 
ranges, while the Morse function is applicable at equilibrium and 
greater separations. In order to better approximate the true 
potential function, various investigators have combined different 
potential functions. As mentioned earlier, EVE [7] used a com- 
posite potential that consisted of an exponentially screened 
Coulomb potential at very close separations, a Born-Mayer po- 
tential at weakly repulsive distances, and a Morse potential over 
the remainder of the potential curve. In their studies, Harrison 
and associates combined a Born-Mayer potential and a Morse 
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potential. These two potentials were joined by a cubic equation 
in the region near their intersection [l,2,3>43. 

C. EXPERIMENTAL RESULTS OF KORNELSEN AND SINHA 

In 1968, Kornelsen and Sinha [ 7 , 8 ] published results concerning 
the binding energies of trapped particles such as neon, argon, 
krypton, and xenon ions in a tungsten surface. The particles were 
forced into the surface from a beam created by an ion gun which 
gave ion energies of 40 eV to 5 keV . The tungsten crystal was 
then heated and the rates of evolution of the trapped gas were 
measur ed. 

The temperature at which desorption peaks occured thus gave 
an indication of the binding energies of point defects in the 
tungsten crystal. Quantitatively sirailiar results were obtained 
with argon, krypton, neon, and xenon ions. Four peaks were ob- 
served below 1650 °K and were labeled as O' peaks. A single peak 
above 1700 °K was measured and was called the p-peak. It was con- 
cluded that the GJ-group of peaks were the results of single step 
desorptions from sites very close to the metal surface. They 
further concluded that the different OJ-peaks could correspond to 
different types of point defect binding energies in the tungsten 
crystal. Specifically mentioned were defects of three types; 

(a) interstitial and substitutional positions in the lattice, (b) 
different distances from the site to the surface, and (c) different 
locations of a nearby lattice defect, such as a vacancy, relative 
to the site and the surface. 
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III. THE SIMULATION MODEL 



A . THE CRYSTAL 

The model used in this research was essentially the same one 
used by Vine Cl4] as explained in Section II. In subsequent dis- 
cussion, when it is necessary to specify certain of the computer 
program variable names, they will be placed in braces. 

All of the simulations in this research were done with tungsten 
crystals of varying sizes. The objective was to use the smallest 
crystal dimensions that would give realistic results. The di- 
mensions described below refer to the number of planes of atoms 
in each of the three rectangular coordinate directions. Simu- 
lations were performed with sizes 8 x 6 x 8 of 96 atoms, 

10 x 6 x 10 of 150 atoms, 10 x 8 x 10 of 200 atoms, and 
10 x 10 x 10 of 250 atoms. Of these, the latter two were judged 
to be of most use because of their greater depth in the y-direction. 
The y direction was always used as the direction of escape for the 
point defect atom. 

Each atom in the simulated crystal was numbered, with the 
first position always assigned to the point defect atom. The 
remainder of the positions were assigned in sequence according to 
the coordinate locations. The numbering was started in the y = 0 
plane and continued until all the atoms in that plane were speci- 
fied. This procedure was repeated for the remainder of the y planes 
in the simulated crystal. 
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B. THE TIMESTEP INTERVAL 



The numerical method of time integration used in the model was 
the average force method [ 13 ]. The value of At {dt} used in this 
procedure was of critical importance in determining whether or not 
the model would approximate reality. Also, the program running 
time was a function of the timestep duration. 

In order to best approximate reality, the {dt} was kept smaller 
early in the program when most movement was expected, and was al- 
lowed to grow larger as equilibrium was approached. The parameter 
that controlled the timestep duration was {Dll}, the distance which 
the most energetic atom was allowed to move in a single timestep. 

In previous work, {DTI} was held constant throughout the duration 
of a program. If the motion was expected to be slow, {dTi} was 
given a larger value than in a situation where simulated motion 
was expected to be greater. For static equilibrium problems, at 
least 100 timesteps were needed to reach an approximately stable 
position. In addition, it became necessary to insert a maximum 
value for the timestep, {dt}, into the program to prevent unrea- 
listic movement of the atoms and breakdown of the model. 

In the current program, changes were made to allow {DTI} to 
decrease during the program. For example, in static runs on He-W, 
Ne-W, Ar-W, and Kr -W, {DTI} was initially given a value of 0.1 
lattice unit {lu} . Tungsten forms a bcc crystal, with LU equal 
to ^LC or 1.58 X, where LC is the Lattice Constant or cube edge 
distance. {Dll} was allowed to change in decrements of 0.01 LU 
per timestep for the first 10 timesteps. Then the {dti} value 
of 0.01 LU was allowed to change in decrements of .001 LU for 



16 



another 10 tiraesteps. Equilibrium was reached by 30 timesteps 
according to this procedure, resulting in a significant decrease 
of computer running time. Also, the equilibrium positions obtained 
in the simulation were closer to the expected (lio) split inter- 
stitial positions than were obtained with a constant { DTl} value. 
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IV. SIMULATION PROCEDURE 



A. STATIC SIMULATIONS 

The first stage of the simulation procedure was concerned with 
finding the equilibrium positions for the point defect of interest 
in the top several layers of a bcc tungsten crystal. Previous work 
done by Johnson [ll] , indicated that a bcc crystal had at least 
two different types of stable split interstitial orientations. 

The first of these was in a (100) direction along a (110) plane. 

The other stable orientation was in a (lio) direction in a (100) 
plane. An analysis of the total of 12 stable split interstitial 
positions possible about a given atom in the two orientations 
listed above, revealed that there were only three independent types 
of sites. (See Figure 1.) The first of these is located on a 
(110) plane with { NVAC} and lies closer to the surface than {nVAc} . 
The second independent position lies in the same (100) plane as 
{NVAC} and is at the same depth in the crystal. A final position 
similar to the first, in that it is along a (110), is located at 
a depth below that of {NVAC}. The three different split inter- 
stitials were postulated to have different binding energies be- 
cause of their varying depths in the crystal. 

In order to reduce the running time in locating the final 
coordinate positions at each interstitial location, the static 
programs were started with the foreign atoms in the approximate 
area of the expected final positions as. discussed above. For 
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Figure 1 



SPLIT INTERSTITIAL... SITES F OR B C C__T UNG STEN 
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example, when it was desired to locate position A, the foreign 
atom was initially placed along the correct (110) plane approxi- 
mately one lattice unit from {nvAC}. 

B. DYNAMIC SIMULATIONS 

The dynamic simulation program used the final positions com- 
puted in the static programs as input initial crystal positions. 
The lattice generator subroutine { B100} , and the point defect 
locator {PLACE} were thus eliminated from the dynamic program. 

Dynamic simulations were broken down into two main categories. 
The first category consisted of survey runs of the areas above 
the interstitials in the direction of escape from the y surface. 
An impact point generator package was inserted into the dynamic 
program, thus permitting the interstitial to be directed at a 
specific number of locations in a predefined area. The results 
of the survey run were analyzed to determine the optimum aiming 
point for the interstitial at a specific initial energy. It was 
observed during the impact testing precedure that the "best" 
aiming points were a function of the initial energy given to the 
interstitial. However, it was also found that the optimum aiming 
points at varying energies were in a generally localized area. 
Thus, once the optimum points were found at the starting energy, 
only a localized region around those points was tested at lower 
initial interstitial energies. The decrementing process of the 
initial interstitial energy constituted the second main phase of 
the dynamic simulations. The procedure in this phase was to de- 
crease the interstitial' energy until it could no longer escape 
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from the crystal* The minimum escape energy was said to be the 
simulated binding energy of the interstitial at the particular 
location tested. 



V. PRESENTATION OF CAT A 



A. STATIC RUNS 

1 • Preliminary Testing of Static Program 

One hundred twelve computer runs were made in the initial 
testing phase of the static program. The bulk of the testing was 
done in four general areas: 1) Program shutdown procedures at 
equilibrium; 2) Optimum crystal size determination; 3) Equilibrium 
position as a function of the initial interstitial position; and 
4) Realistic timestep determination procedure. 

The best method of stopping the program at equilibrium has 
been a matter of concern for many years with the static simulation 
program. In the present investigation, the first method tested was 
a shutdown procedure initiated whenever a sharp { DT} decrease was 
encountered in the program. This test proved to be of limited 
success, and was later abandoned. Another method that was at- 
tempted was a test of {emax} against a value such as . 04 , to deter- 
mine if equilibrium had been reached. This method also proved only 
partially successful. Finally, a test was made of the average 
kinetic energy of an atom in the simulated crystal. The average 
energy was taken to be the total kinetic energy [tPKe} divided by 
the number of atoms (ll} , or equivalently, {TPKE} multiplied by 
the reciprocal of the number of atoms {rll}. The crystal was 
assumed to be at equilibrium if the average kinetic energy of an 
atom was less than or equal to the value 0.025 eV, a value for the 
average thermal energy of an atom. Satisfactory results were 
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sometimes obtained by this method but the test had to be removed 
in many cases to allow the simulated crystal to run a greater 
number of timesteps and reach equilibrium. 

As previously mentioned, different sizes of crystals were 
simulated to determine the optimum dimensions for the crystallite. 
Many tests were made on crystal sizes smaller than the 10 x 10 x 10 
used by Vine Cl4]« Although smaller sizes such as 10 x 6 x 10 often 
gave reliable results, it was finally decided to use the 10 x 10 x 10 
size for the reported split interstitial positions in order to have 
a standard size applicable over all crystal positions of interest. 

In his work, Vine placed his interstitials in the middle 
of the open channels in the simulated crystal. Numerous runs in 
the present research have indicated that equilibrium is reached 
more quickly when the initial starting positions are not in 
channel centers, but along the directions of the expected splits 
in the general area of the final positions. 

All of the initial work was done with a fixed {DTI} value 
in the program. The [DTl} value normally used was .05 LU. The 
timestep interval was later modified as explained in The Simulation 
Model and the computer running time was cut by approximately two- 
thirds due to this procedure. 

2 . Positions for Helium in Tungsten 

Helium was the lightest point defect used in the static 
simulation model. It was thought that the small, light atom would 
essentially do all the movement and come to rest in a channel 
center \J~2~ LU away from {nVAC} along the (lio) direction. The 
results of the runs show that the movement was almost as expected. 
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Also, a comparison of the C site for atom 64 and the A site for 
atom 114 indicates that there is a strong possibility that these 
two sites are degenerate. This is based on the fact that the 
potential energies of both sites are close, and also that the 
[nVAc] for each site is displaced only a negligible distance. The 
same probability of degeneracy is also seen to exist for C -89 and 
A-139. The determination of degeneracy of sites is seen to be a 
complex evaluation of potential energy differences, movements of 
{nVAC} , and the potential gradient between the two sites. The 
final positions for the A, B, and C split interstitial positions 
in the third, fourth, fifth, and sixth layers where [NVAC] was 
64, 89 , 114, and 139 respectively are shown in Figures 2, 3 and 4- 
The results obtained with helium were thought to be somewhat in 
error, since the heavier atom, neon, moved further in its simu- 
lations than helium. Further work is needed to confirm the equi- 
librium positions for helium. 

3 . Positions for Neon in Tungsten 

The neon split interstitial locations were simulated for 
the same positions as helium. (See Figures 5, 6 and 7.) Essenti- 
ally, the {NVAC} atom remained in its lattice site and the neon 
moved along a (lio) direction to a distance \f!T LU away from 
1 NVAC} . 

4 • Positions for Argon in Tungsten 

The bulk of the testing of this research was done with 
argon as the simulated point defect. The results of the runs for 
the 10 x 10 x 10 crystal size are shown* in Figures 8, 9 and 10 . 
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Figure 4 _ 

B SITES FOR ATOMS 64, 8 9, 114, 1-39 
HELIUM — Q TUNGSTEN— Q) 
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Figure 5 
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Figure 7 

B SITES FOR ATOMS 64,89, 114,13a 
NEON-0 . TUNGSTEN -Q) 
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Figure 8 
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Figure 9 
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Figure 10 



B SITES FOR ATOMS 64, 89,114,139 
ARGON -Q TUNGSTEN -Q 
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The expected (lio) splits were observed with {NVAC} moving away 
from its site in relation to the f inal pos it ion of the argon defect. 



5 • Positions for Krypton in Tungsten 

The split interstitial positions for krypton in a simu- 
lated tungsten crystal were also plotted. (See Figures 11, 12 
and 13*) Once again, satisfactory results were obtained and the 
(lio) splitting was observed. 

6 . Positions for Xenon in Tungsten 

For the Xenon Runs, the tungsten repulsive potential 
function was used for the xenon potential. The initial runs, 
using the {dti} employed for the other defects, failed to give 
the postulated split interstitial positions. At least two factors 
were seen to complicate the Xenon-Tungsten runs. Firstly, the 
mass of Xenon was approximately 6/7 that of tungsten, so the 
lattice site would be shared almost evenly. This would mean that 
tungsten would have to move a significant distance from its lat- 
tice position. Secondly, the relatively large size of the xenon 
defect would require more movement of the surrounding atoms to 
accommodate the extra atom. 

The {dti} scheme was changed to allow for more movement of 
the atoms by allowing {dti} to change in smaller decrements. 
Another method that was used was an initial displacement of both 
the xenon and the tungsten away from the lattice site. When both 
atoms were displaced, it became necessary to use an initial {dti} 
of .05 LU or less. An average value of the C site for atom 89 
was computed. (See Figure 14 . ) More work is needed to simulate 
all of the interstitial positions. 
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Figure 12 



A a C SITES FOR ATOMS 89 , 139 
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Figure 13 

B SITES FOR ATOMS 64, 89, 114, 139 
KRYPTON -O ... TUNGSTEN -O 
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7 . Split Interstitial Distance Ratios as a Function of 



Relative Masses 

During the static testing phase, the ratios of the split 
interstitial distances from the initial {iWAc} were measured. An 
attempt was made to correlate these ratios to the inverse ratios 
of the atomic masses of the atoms involved. (See Table I.) The 
point defect results of helium and neon did not give any signi- 
ficant correlation, but the Argon and krypton atoms gave ratios 
of splits in good agreement with the expected values from mass 
ratio calculations. 

B. DYNAMIC RUNS 

1 . Survey of the Possible Directions of Escape for the Ar -W 

Simulated Split Interstitial in Site A-8Q 

As was mentioned in the Section IV, survey runs were made 
of the possible escape directions in the plane one LU. above the 
defect. The results for the survey for the Ar -W split interstitial 
in site A -89 were plotted. (See Figure 15.) [ Dy} , the distance in 

lattice units that the argon moved in 25 timesteps, and {vy}, the 
velocity that the argon had after 25 timesteps are shown for each 
impact point. Due to the multiple scattering of the defect as it 
moved toward the surface, 25 timesteps were not sufficient to 
provide conclusive evidence at any one impact point. The survey 
was limited to 25 timesteps per impact point, however, due to 
computer time considerations. The main benefit gained from the 
survey run was the elimination of certain areas from further testing, 
such as those on the outer perimeter of the survey area. Also, much 
information was gained on the most likely area of escape for the 
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TABLE I 








Theoretical Splits 
From Lattice Site 


Simulated Splits From 
La tt ice Site 




^Tungsten Mass A 
\ Impurity Mass / 


/ Impurity Distance \ 
(Tungsten Distance j 


He-W 


45.96 




Ne-W 


9.11 




> 

K 

1 

£ 


4.60 


4.36 


7 * 

K 

1 

£ 


2.19 


2.18 


Xe-W 


1.40 
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Figure 15 



IMPACT £REA : DYNAMIC SIMULATION OF ESCAPE 
DIRECTIONS FOR ARGON IN SITE A- 89 



/ DY=-0.31 \ 


/ dy=-i.73 \ 


/ DY= -1.28 \ 


/dy=-i.73 \ 


/ DY= -1 . 32 \ 


5.99 [ vy=10500 L 


vy=1954 u 


_/ vy =6050 L 


_/ VY=-3726 L 


-1 VY=1330 ] 



/ dy=-0.37 \ 


/dy=-1.33\ 


f DY= - 1 . 2 2 \ 


/dy=-1.50\ 


/ dy=-o.85\ 


5.5 9 ( VY=7970 L 


J VY=7844 L 


J VY=8937 ]_ 


_J VY=5000 L 


J VY=383 ) 



5.19 



’dy=-o.i4< 

vy=8476 



dy=-i .48 

VY=3003 



DY=-1 . 33 
vy=7510 



dy=- 0.98 

VY=-l420 



dy=o.44 

vy=8ii7 



dy=-o. 13 '; 
4.79 [ vy=8672 



DY=-1.28/ 

vy=3916 



dy=-i.34 

vy=7551 



dy=-o.95 

vy=1355 



DY= . 4 50 
VY=8091 



/ dy=-o.34 \ 


/ DY=-1.37 \ 


/ DY = - 1 .23 \ 


/ dy=-i.47\ 


/ dy=- 0.67 \ 


4.39J vy=8578 J- 


-I vy=7738 }- 


_{ vy=8678 J- 


J vy=3795 


-I VY=2159 ) 



DY=-0.30r 
3.99[ vy=10646 



3.44 



\ / DYZ 


: -i . 71 \ 


/ DY= - 1 . 2 9 \ 


/ DY=-1 . 53 \ 


J— ( VY= 


: 2i45 U 


J VY=5907 V- 


_| vy= 657 l 


3 


'“84 


47?4 


4T64 



DY=-1. 13 
VY=279 5 



5 .04 



41 



defect. From the results of the survey run, and from symmetry 
considerations of the open channel, it was decided to concentrate 
further testing on several points with the same Z coordinate as 
the interstitial atom. 

2 . Determination of the Simulated Binding Energy of Argon 
in Site A-89 

Detailed testing of several impact points with the same Z 
coordinate as the argon was carried out for 100 timesteps per 
impact point. This was generally enough time for' the simulated 
interstitial to escape the crystal it its initial energy was 
great enough. The results are summarized below. 

TABLE II Detailed Impact Point Testing 



Impact Points in 
Plane 1.0 LU in 
-y direction from 
Ar gon 


Initial Energy 


20 eV 


10 eV 


8 eV 


6 eV 


4 eV 


Cx = 3.34 LU 
CZ = 4.99 


Yes 


Yes 


No 


No 


No 


cx = 4.24 
CZ = 4.99 


Yes 


Yes 


Yes 


Yes 


No 


cx = 4.44 

CZ = 4.99 


Yes 


Yes 


Yes 


Yes 


No 



Yes - Atom Escapes Crystal 
No - Atom Does Not Escape 
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VI . CONCLUSIONS 



The method used to decrease {dti} on every timestep appears 
to be a significant improvement over the old method of keeping 
{Dll} fixed. In effect, the new method forces the simulated 
crystallite to come to equilibrium in a predefined number of 
timesteps. In addition to a saving of computer time, this method 
also eliminates the need for a equilibrium shut down procedure 
in the program. The same {dti} decrement scheme could not be used 
for all the point defects tested in this research. For example, 
Argon and Krypton were able to come to their expected equilibrium 
positions using the scheme described in Section III, but Xenon 
and Helium were not. Thus it appears that defect size, as well 
as the expected degree of movement may dictate the [dti] decre- 
ment scheme. 

The expected (lio) split interstitials were confirmed for the 
helium, neon, argon, krypton, and xenon foreign defects. The re- 
lative degree of splitting does appear to be related to the masses 
of the interacting atoms, but conclusive evidence of this was only 
obtained for the Argon-Tungsten, and Krypton-Tungsten runs. 

The dynamic runs have shown the direction of escape from the 
crystal to be a mult i-coll is ional process, with the defect under- 
going many intermediate direction changes at the low kinetic 
energies tested. The most likely escape direction was also seen 
to vary somewhat with the initial kinetic energy given to the 



defect . 



The simulated value of approximately 4 eV for the binding 
energy of Argon in site A-89 appears to be a reasonable value. 
Kornelsen ! s data showed an energy level at this value and his 
levels were postulated to arise from defects in the first several 
layers of the tungsten crystal. 

Further work is needed to confirm the remainder of the energy 
levels found by Kornelsen, by testing other point defect sites in 
the top layers of the simulated crystal. 
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APPENDIX A: COMPUTER PROGRAM GLOSSARY 



NOTE: In this glossary, the terms "point defect atom", "bullet", 

and "primary" are synonymous; and the terms "latttice atom" and 
"target" are synonymous. 

AC: Distance measurement used in impact point generator 

ALPHA: Input Morse potential parameter 

BSAVE: Target mass/(target mass + bullet mass); distributed 

potential energy between target and bullet 
BIND: Negative of the total potential energy (TPOT) at time 

zer o 

EMAS: Mass of bullet in amu 

BULLET: Alpha-numeric array for point defect material 

CFO, CFl , CF2 : Force parameters of cubic fit between Morse and 

Bor n -Ma yer functions 
CGB1, CGB2: Morse potential parameters 

CGD1 , CGD2: Morse potential parameters 

GGFl , CGF2 : Morse force parameters 

COX, COY, COZ: Cosines of angles to x,y, and z axes respectively 

CPO, CPI, CP2 , CP3 : Potential parameters of cubic fit between 

Morse and Born-Mayer functions 

CVD: CVR x 10 ^ , converts lattice units to meters 
, -19 

CVE: 1.6 x 10 , converts electron volts to joules 

CVED: CVE/CVD, a ratio used to avoid repeated division 
-27 

CVM: 1.672 x 10 , converts atomic mass units to kilograms 

CVR: LU in angstroms; converts lattice units to angstrom units 
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CX, CY, CZ: Coordinates of impact point 

D1X, D1Y, D1Z: Displacement coordinates for location of inter - 





stitial from reference atom, NVAC 


DCON : 


Input Morse potential parameter 


DDT I: 


Time increment that is subtracted from DTI after each 
t imes tep 


DFF : 


ROE-DIST, the distance closer than ROE that an atom is to 
the primary 


DIST: 


Distance between any two atoms 


DLPE: 


TLPE-TLPE0, the change in total local potential energy 




since time zero 



DRX, DRY, DRZ : x,y,z components of DIST 



DT: 


Length of a timestep in seconds 


DTI: 


Number of lattice units most energetic atom may move in 
one timestep 


DTIS: 


Starting value of DTI 


DTOD: 


DT/CVD--a ratio used to avoid repeated division 


DTOM: 


DT/PTMAS--a ratio used to avoid repeated division 


DTOMB: 


DT/PEMAS--a ratio used to avoid repeated division 



DX(I), DY(I), DZ(I): Change in position of ith atom from initial 





position at time zero 


EMAX: 


The maximum energy encountered in any cycle 


ERAT: 


Measure of the average kinetic energy of an atom 


EV: 


Primary energy in electron volts 


EVR : 


Primary energy in kilo-electron volts 



EXA, EXB: Input Born-Mayer potential function parameters for the 

target 
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F2: 


Square of the force on a specific atom 


FA: 


The component force increment on an atom 



FDTI: DTI X CVD, a parameter used to determine DT by maximum 





energy method 


FM: 


A small number used in checking potential energy zero 
point 


FM2: 


FM squared 


FMAX: 


Maximum total force on the most stressed atom in the 
crystal 


FMAX1 : 


Maximum total force on Atom 1 


FOD: 


FORCE/DIST--a ratio used to avoid repeated division 


FORCE : 


Numerical value of the force function with a variable 




parameter 



FX(I), FY(I), FZ(I): x,y,z components of total foce on an atom 



FXA: 


Born-Mayer force function parameter 


HBMAS: 


% BMAS-- a ratio used to avoid repeated division 


HDTOD: 


\ DTOD-- a ratio used to avoid repeated division 


HDTOM : 


\ DTOM-- a ratio used to avoid repeated division 


HDTOMB: 


\ DTOMB-- a ratio used to avoid repeated division 


HTMAS : 


\ TMAS-- a ratio used to avoid repeated division 


11 : 


Variable in cubic fit subroutine 


13 : 


Variable in cubic fit subroutine 


I DEEP: 


Number of mobile layers 


IHl : 


Alpha numeric array for program title 


IH2 : 


Alpha numeric array for Morse function parameters 


IHB: 


Alpha numeric array for bullet element 


IHS: 


Alpha numeric array for type and orientation of crystal 
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IHT : 


Alpha numeric array for target element 


ILAY : 


Same as I DEEP 


IN: 


Odd-even integer used to determine atom site establishment 


IP: 


Subscript value of atom. Used in subroutine STEP and 
ENERGY 


IQ: 


Parameter that determines whether or not a self defect is 
to be given a repulsive potential or a composite attractive- 
repulsive potential 


I SHUT: 


A parameter used to shut down the program 


IT: 


Unsealed fixed point x coordinate used in lattice generation 


ITT: 


Odd-even integer used to determine atom site establishment 


ITYPE: 


Parameter used to determine the type of point defect: 
vacancy, self-interstitial, replacement, foreign inter- 
stitial 


IVACX, 


IVACY, IVACZ: Input plane numbers to specify NVAC 



IX, IY, IZ: Number of x,y,z planes of crystal 



J2 : 


Variable in the cubic fit subroutine 


JT: 


Unsealed y coordinate used in crystal generation 


JTS : 


Variable used to establish atom sites 


JTT : 


Variable used to establish atom sites 


KF: 


Final K in LOCAT (K) assigned to an atom 


KT: 


Unsealed z coordinate used to establish atom site 



LCUT(I): Used to identify an i th atom which is not included in 





calculations 


LD: 


The highest numbered atom in the mobile layers 


LL: 


The highest numbered atom in the entire crystal 
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LOCAT(K) : Dimensioned variable that remembers the numbers of the 





atoms within a radius ROEL of the primary at time zero 


LS: 


Variable associated with each of the nine lattice gene- 
rator subroutines 


MCRO: 


One number higher than the order of the fit between the 
Born-Mayer and Morse potentials, always 4 in this simu- 
lation 


ND: 


Data output increment, in numbers of times teps 


NDEC : 


Counting index for DTI variation 


NEW: 


Parameter used to determine whether or not atom numbers 
have been stored in LOCAT(K) 


NPAGE : 


Page numbering variable 


NR UN: 


Parameter used to determine whether or not to read 
additional data cards 


NS: 


Initial print statement tiraestep number 


NT: 


Timestep number 


NTT: 


Timestep number limit before shutdown 


NVAC: 


An atom number used to establish point defects or used as 
a reference point for interstitial placement 


PAC: 


Parameter for bullet force function correction 


PEMAS: 


Primary mass in kilograms 



PEXA, PEXB: Input Born-Mayer potential function parameters for the 





bullet -tar get interaction 


PFPTC: 


Primary force function evaluated at ROE 


PFXA : 


Primary force function parameter 


PKE( I) : 


Kinetic energy of the ith atom 


PLANE: 


Alpha-numeric array for lattice orientation 
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POT: 


Potential energy between two atoms 


PPE( I ) : 


Potential energy of the it_h atom 


PPTC: 


Primary potential function evaluated at ROE 


PTE( I) : 


Total energy of the itti atom (potential + kinetic) 


PTMAS : 


Target mass in kilograms 


RE: 


Input Morse potential parameter 


RLL: 


Reciprocal of LL 


RO: 


Spacing constant in FCC(llO) lattice generation subroutine 


ROE: 


Nearest neighbor distance 


ROE2 : 


ROE squared 


ROEA: 


Maximum cut off for Born-Mayer potential 


ROEB: 


Minimum cut off for Morse potential 


ROEC: 


Maximum cut off for Morse potential 


ROEC2 : 


ROEC squared 


ROEL: 


Radius inside of which local potential energy is found 


R0EL2 : 


ROEL squared 


ROEM: 


ROE-DTI, region in which modification of repulsive force 




must be made 



RX(I), RY ( I ) , RZ(I ): x,y,z coordinates of an ith^ atom at any time 
RXI(I), RYI(I), RZI(I): x,y,z coordinates of an ith_ atom's initial 
position 

RXK(I), RYK(I), RZK(I): x,y,z coordinates of temporary position of 





an iLh atom during force cycle 


SAVE: 


^ POT 



SCX, SCY, SCZ: x,y,z coordinate scale factors 



START : 


An optional timing variable, not used in this simulation 


SUM: 


Variable in cubic fit subroutine 
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TARGET : 


Alpha -numer ic array for target material 


TSAVE : 


Bullet raass/( tar get mass + bullet mass); distributes 



potential energy between target and bullet 



TE: 


Total energy of all crystal atoms (kinetic + potential) 


TEMP: 


Temperature of lattice in degrees Kelvin. Not used in 



this simulation 



TEFAC: 


Product of DTI * CVD 


TFAC: 


A time factor ratio used to determine DT by maximum force 



method 



TFACB: 


TFAC for the bullet 


THERM : 


Thermal energy of atom. Not used in this simulation 


TIME: 


Elapsed problem time in seconds 


TLPE: 


Total local potential energy of atoms within a radius ROEL 


TLPE0 : 


TLPE at time zero 


TMAS: 


Target atom mass in amu 


TPKE: 


Total kinetic energy of all crystal atoms 


TPOT : 


Total potential energy of all crystal atoms 


TPOTL : 


Storage position for the last computed value TPOT 


VSS : 


Storage variable for velocity components 



VX(I), VY(I), VZ(I): x,y,z components of itji^ atoms velocity 



X, Y, Z: 


Unsealed coordinates used in crystal generation 


XSTART : 


X Coordinate used in impact point generator 


YLAX( I ) : 


Relaxation in -y direction of ith_ layer in L.U. 


ZP: 


Floating point form of JTT 
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COMPUTER PROGRAM 



THE FOLLOWING PROGRAM IS THE ONE USED FOR THE STATIC 
SIMULATION RUNS. THE DYNAMIC PROGRAM IS ESSENTIALLY THE 
SAME AS THE STATIC. IN THE DYNAMIC PROGRAM, SUBROUTINES 
CBIOOJ AND (PLACE) ARE OMITTED, AND THE INITIAL ATOM POSI- 
TIONS ARE READ IN ON DATA CARDS. ALSO, AN IMPACT POINT 
GENERATOR PACKAGE IS INCLUDED IN THE DYNAMIC PROGRAM. 

BOTH PROGRAMS CALCULATE THE EQUATIONS OF MOTION FOR THE 
SYSTEM THROUGH THE USE OF THE APPROPRIATE POTENTIAL FUNCTION 
AND THE AVERAGE FORCE METHOD OF INTEGRATION. THE NET 
DISTANCE OF MOVEMENT, VELOCITY, AND ENERGY VALUES ARE 
PRINTED OUT FOR THE DESIRED ATOMS ON SELECTED TIMESTEP 
NUMBERS. THE FINAL POSITION, VELOCITY, AND ENERGY VALUES 
FOR EACH ATOM ARE SUMMARIZED AT THE END OF THE PROGRAM. 

// EXEC FORTHCLG, TIME. G0=20, REGION. GC=140K 
//FORT . SYS I N DD * 

C 

DIMENSION VX(IOOO) ,VY( 1000) ,VZ( 1000) ,PKE( 1000) 
DIMENSIONING OF VARIABLES NOT NEEDED IN COMMON 
C 

DIMENSION DX (1 000) ,DY( 1000) ,DZ( 1000) ,PTE( 1000) 
DIMENSION RXKdOOO ) ,RYK{ 1000) , RZK( 1000 ) 

COMMON LABELING OF VARIABLES REQUIRED IN OTHER SUBROUTINES 
COMMON/ COM l/RXl 1000) ,RY( 1 000 ) , R Z ( 1 000 ) ,LCUT( 1000) , 

ILL ,LD, I TYPE , NVAC 

COMMON /COM 2/ 1 H 1 ( 20 ) , I H2 ( 8 ) , I H S ( 1 0) , I HB ( 6 ) , I HT ( 6 ) , 
1TARGET (4) , TM AS , BULL ET ( 4 ) , BM AS , P LAN E , TE MP , THE RM 
COMMON /COM 3/ RX I ( 1000) , RY I ( 1 000 ) , RZ I ( 1000 ) , C V R , EVR, 

1NT, TIME, DT, DTI , I LAY 
1 IVACX, IVACY, IVACZ 

COMMON /C0M4/ I X,IY,IZ,SCX,SCY,SCZ,IDEEP,D1X,01Y,D1Z, 
COMM ON /COM 5 /RO E , RO C 2 , R OEM , E XA , E XB , PE XA , P EXB , F XA ,P F X A , 

1 IQ,TSAVE, BSAVE 

COMMON /C0M6/FX ( 1000) ,FY(1000) ,FZ(1000) ,PAC,PFPTC, FM 
COMMON/ C0M7/P PTC ,T POT, PPE( 1000) , TLP E ,R OEL , R0EL2 ,NEW 
C0MM0N/C0M8/R0EA, ROEB,ROEC , RO EC2 , C PO , C PI , CP 2 , CP3 , 

1CF0, CF 1,CF2,CGD1 ,CGD2,CG61 ,CGB2,CGF1 »CGF2 
COMMON/COMA/ A(4,5),MCR0 
C 

READ STATEMENT FORMATS 

9010 FORMAT (20A4) 

9020 FORMAT ( 8A4,5F8. 5»2F5. 2) 

9030 FORMAT ( 4A4 , 3 F8 . 5 , 6 A4 , F 6. 2) 

9040 FORMAT ( Fo.3,5X, I5,6I4,3F5,3 12 ) 

90 50 FORMAT ( 10A4, A4 , 4 1 3 , F 8. 4 , 1 4 , F 5 . 3 ) 

C 

WRITE STATEMENT FORMATS 
9610 FORMAT ( 1 H 1 ) 

9620 F0RMAT(47X, 'SUMMARY OF ATOMS ' / / , 35 X , 8 A4 , • , NT ='I4,//, 
1 3 ( • ATOM POSITION BIND ENERGY •),//) 

9630 FORMAT (3 ( 15,31=6.2, F3. 4, 8X) ) 

9640 FORMAT (/4X,F10. 3, 25H EV, TOTAL KINETIC EN ERGY , , F 10 . 3 , 

1H EV, TOTAL POTENTIAL ENERGY ,F 10.3 , « EV, REDUCTION', 
1//60X' RADIUS =' , F5 .2, ) 

96 50 FORMAT ( 1 05X , 4H P AGE , I 3 , / , 1H 1 ) 

9660 FORMAT ( / • ATOM DX DY DZ 

1VX VY VZ KE PE TE'/) 

9670 FORMAT ( 1 1 3 ,3F 1 0. 3 , 3F 10. 1 ,3F10.4 ) 

9680 FORMAT ( ' SHARP DT DECREA S E ’ , 2E 1 0 . 3 ) 

9690 FORMAT (I4,3F5.2,I4) 

9691 FORMAT ( 9F8 . 4 ) 

9692 FORMAT ( IX, 14,/ ) 

INITIALIZING 

START=0.01*ITIME(XX) 

DO 2 1=1,1000 
RXK ( I ) =0 . 0 
RYK ( I ) =0.0 
RZK ( I ) =0. 0 
VX ( I ) = 0 . 0 
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2 



VYl I )=0.0 
VZ( I ) = 0.0 
PKE ( I ) =0.0 
PP E ( I ) =0 .0 
PTE ( IJ =0.0 

rzi ( n=o.o 

I SHUT = 1 
NRUN=0 

c 

INPUT DATA 

READ ( 5,9010) 

READ ( 5,9020) 

READ ( 5,9030) 

READ ( 5,9030) 

READ ( 5,9050) 

C 

CONSTANTS AND SCALING 

DTI S=DT 1 
RCE2=4 . 0 
ROE=S QRT ( ROE 2 ) 

ROEM = ROE-DTI 
RUEL2 = R0EL S ' : R0E L 
CV E= 1 . 60 E- 19 
CVM=1.672E-27 
VF AC = . 50 
FM= 1 . 0 E-10 
F M2 = F M* F M 
CVD=CVR*1.0E— 10 
CVED=CVE/CVD 
PTMAS=TMAS*CVM 
PBMA S= 8MA S~C VM 
HT MAS = 0 . 5*PTM AS / C.V E 
HBMAS=0. 5*PBMAS/CVE 
TSAVE=BMAS/( BMAS+TMAS) 

BSAV E=TMAS/ ( BMAS+TMAS) 

C 

REPULSIVE POTENTIAL PARAMETERS 

FXA=ALOG(— EXB*CVED)+EXA 
PFXA=ALOG(-PEXB*CVED)+PEXA 
PPTC=EXP (PEXA+PEX3*R0E ) 

PAC=ALOG(CVED) +PEXA 
PF PTC= E XP ( PAC+PEXB*ROE ) 

C 

ATTRACTIVE POTENTIAL PARAMETERS 

CGD1=AL0G( DCONJ+2. 0*ALPHA*RE 
CGD2=AL0G(2. 0*DCON)+ ALPHA* RE 
CGBl=-2. 0*ALPHA*CVR 
CGB2=— ALPHAS CV R 
CGF 1 =A LOG I -CGB 1 *C V ED ) +CGD1 
CGF2=AL0G(-CGB2*CVED) +CGD2 
C 

CUTOFF DISTANCES FOR ATTRACTIVE AND REPULSIVE POTENTIALS 
ROE A= 1 . 50/CVR 
ROE 8=2 . 0/ CVR 
R0EC2=R0EC*R0EC 

C 

PARAMETERS FOR CALCULATION OF THE BEST CUBIC FIT IN THE GAP 
BETWEEN MAXIMUM DISTANCE CUT O p F OF THE REPULSIVE POTENTIAL 
(ROEA), AND MINIMUM DISTANCE CUTOFF OF THE ATTRACTIVE POTEN- 
TIAL ( ROEB ) . SUBROUTINE CROSYM ACTUALLY PERFORMS THIS CURVE 
FITTING. 

A(l,l ) = 1 .0 

A { 1 , 2) =ROE A 

A(1,3)=R0EA*R0EA 

A(1 ,4) =R0EA**3 

All, 5) =EXP(EXA+EXB*ROEA) 

A( 2, 1 ) = 1.0 



I H 1 

IH2,DC0N, ALPHA, RE, ROEC,ROEL 
BULLET , BMAS , PEXA ,P EXB, I H8, THERM 
TARGET, TMA S, E XA,EXB, I HT, TEMP 
IHS,PLANE,LS,IX, 1Y, IZ,CVR,MCRC ,DTI 

FACTORS 
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A(2 » 2 ) = ROEB 
A( 2,3) =ROEB*ROEB 
A( 2» 4) =R0EB**3 

A ( 2 , 5 ) =EXP (CGD1+CG31*R0EB)-EXP(CGD2+CGB2*R0EB) 

A( 3, 1) =0.0 

A ( 3 , 2 ) =— 1.0 

A(3,3)=-2.0*R0EA 

A( 3, 4) =-3 . 0*ROEA*ROEA 

A ( 3 , 5)=EXP(FXA+EX8*R0EA)/CVED 

A { 4 , 1 ) =0.0 

A ( 4 , 2 ) =-1.0 

A ( 4 , 3 ) =— 2 . 0*R0 E 3 

A (4, 4) =-3.0*R0EB*ROEB 

A( 4, 5) =(EXP(CGF1 + CGB1*R0EB)-EXP ( CG F2 + CGB2 * ROE B ) J/CVED 
CALL CROSYM 
CPO=A ( 1 , 5 ) 

CP1=A( 2,5) 

CP2=A( 3,5) 

CP3=A (4,5) 

CF0=-CP1*CVED 
CF1 =-2 . 0*CP2*CVED 
CF2=-3 . 0*CP3-CVED 

5 READ ( 5,9040) EVR , NTT , NS , ND , I P , I DE EP , I TY PE , NVAC , D IX 
1D1Y,D1Z, IVACX, IV AC Y , IVACZ 
IF(NTT.EQ.O) GO TO 9999 
I Q= I TYPE-1 
EV=EVR*1 .OE+3 

SELECTION OF THE DESIRED CRYSTAL STRUCTURE AND ORIENTATION. 
100, 110, AND 111 PLANES OF FACE-CENTERED, BODY-CENTERED, 
AND DIAMOND STRUCTURES ARE ALLOWED. I LAY AND IDEEP ARE VAR 
IABLES ESTABLISHING THE NUMBER OF MOBILE LAYERS IN THE 
CRYSTAL. RXI(I) AND P.XK ( I ) ARE VARIABLES SAVING THE ORIGIN 
AL X- POSITION OF THE I ' TH ATOM. Y AND Z POSITIONS ARE 
ANALOGOUS. 

14 CALL B100 
30 I LAY= I DEEP 

IF(IDEEP) 35,35,40 
35 LD=LL 
ILAY=I Y 
40 RL L = 1 . O/LL 
TP0TL = 1 . 0 
DO 45 1=1, LL 
RXK ( I ) =RX ( I ) 

RYK ( I ) =RY( I ) 

RZK ( I ) =RZ ( I ) 

RXK I ) =RX ( I ) 

RYI ( I ) =RY( I) 

45 RZI ( I ) =RZ ( I ) 

C 

THIS SECTION ALLOWS ONE TO REPEAT A RUN OF THE PROGRAM WITH 
DIFFERENT DATA WITHOUT REPEATING INITIALIZATION, POTENTIAL 
PARAMETER CALCULATIONS AND CRYSTAL LATTICE BUILDING. SUB- 
ROUTINE PLACE USES LCUT(I) AND NVAC TO CREATE VACANCIES, 
INTERSTITIALS, AND REPLACEMENT IMPURITIES AT DESIRED LOCA- 
TIONS IN THE LATTICE. 

IF(NRUN.EQ.O) GO TO 60 
DO 55 1=1, LL 
LCUT ( I ) =0 
RX ( I ) = RX I ( I ) 

RY ( I ) =RY I ( I ) 

RZ ( I ) =RZ I ( 1 5 
RXK ( I ) =RXI (I) 

RYK ( I) =RYI I I ) 

55 RZK ( I ) =RZ I (I ) 

60 NRUN= 1 

CALL PLACE 
RXK 1) =RX( 1) 

RYI ( 1 ) =RY( 1) 

RZI ( 1) =RZ( 1) 
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RXK ( 1) = RX( 1) 

RYK( 1 ) =RY ( 1) 

RZK ( U =RZ( 1) 

DO 65 1=1, LL 
VX ( I )=0.0 
VY( I J =0. 0 
VZ ( I )=0.0 
PPE ( I )=0.0 
PKE ( I) =0.0 
PT E ( I ) =0. 0 
TPOT =0 . 0 
NE W = 0 

c 

THE ENERGY SUBROUTINE CALCULATES THE POTENTIAL ENERGY OF 
EACH ATOM IN THE LATTICE. SUBROUTINE LOCAL SUMS UP THIS 
ENERGY FOR ALL ATOMS WITHIN A SPECIFIC RADIUS OF THE POINT 
DEFECT. 

THIS SECTION PRINTS OUT X, Y, AND Z COORDINATES, IN LATTICE 
UNITS, AND BINDING ENERGIES OP EACH ATOM IN THE CRYSTAL AT 
TIME ZERO. 

CALL ENERGY 
B I ND=-TPOT 
TE=TPCT+8 I ND 
C 

TIME=0.0 
NT =0 

WRITE ( 6,9610) 

WRITE ( 6,9620) IH2,NT 
DO 70 1=1, LL, 3 
K= I + 1 
J= I +2 

70 WRITE ( 6,9630) I , RX ( I ) , RY ( I ) , R Z( I ) , P PE ( I ) , K , RX ( K ) , 
1RY ( K ) , RZ ( K ) , PPE IK) ,J,RX( J) ,RY(J),RZ(J) ,PPE(J1 
WRITE ( 6 , 9640 ) TPKE, T-PJDT , TE , ROEL 
NPAGE=1 
NPAGE=NP AGE+ 1 
WRITE ( 6,9650) NP AGE 
C 

THIS IS THE MAIN BODY CF THE PROGRAM. BY USE OF THE AVERAGE 
FORCE METHOD, EXPLAINED IN REF. 13, IT DOES ALL 
THE DYNAMICS FOR EACH INDIVIDUAL ATOM. SUBROUTINE STEP 
CALCULATES ALL MUTUAL FORCES AMONG THE ATOMS. BASED ON THE 
FORCES, THIS SECTION THEN CALCULATES TEMPORARY POSITIONS FOR 
THE PRIMARY, AND ALL OTHER ATOMS; RECALCULATES FORCES IN 
STEP; AND THEN RECALCULATES FINAL POSITIONS FOR THE PRIMARY 
AND ALL OTHER ATOMS, BASED ON THE AVERAGE CF THESE TWO 
FORCES. THIS SECTION ALSO INCLUDES ALL KINETIC ENERGY CAL- 
CULATIONS, BASED ON THE VELOCITIES INVOLVED; AND FINALLY 
CALCULATES A NEW TIMESTEP DURATION FOR USE IN THE NEXT TIME- 
STEP, BASED ON EITHER A MAXIMUM ALLOWED FORCE, OR MAXIMUM 
ALLOWED ENERGY. VFLGCITIES ARE HALVED AT THE END OF EACH 
TIMESTEP AS A METHOD OF DAMPING. 

DT = 1 . 0 E— 15 
DDT I =. 005 
NDEC=0 

100 DTI=DT I-DDTI 
NDEC =NDEC+ 1 
DTOD=DT/CVD 

TFAC=2 .0*PTMAS*DTI*CVD 
TFACB=2.0*PBMAS*DTI*CVD 
TE FAC = DT I *CVD 
HDT0D=0.5*DTCD 
DTOM=DT /PTMAS 
HDT0M=0.5*DTC'M 
DT 0M3= DT/P BMAS 
HDTOMB = 0. 5*DT0MB 
200 CALL STEP 

IFILCUT (1) .GT.O) GO TO 240 
1 = 1 

RXK ( I ) =RX ( I ) 



55 



RYK ( I } =RY ( I ) 

RZK ( I ) =RZ ( I ) 

RX ( I ) =RX ( I )+DTQD* { HDT0M8 : -FX ( I ) + VX { I ) ) 

RY ( I ) = RY ( I )+DTGD*( HDT0M8*FY( I )+VY( I ) ) 

RZ( I ) =RZ( I )+DTOD*(HDTOMB*FZ( I ) +VZ( I ) ) 

240 DO 245 1=2, LD 

IF(LCUT ( I ) .GT .0 )G0 TO 245 
RXK ( I ) =RX ( I ) 

RYK ( I ) =R Y ( I ) 

RZK ( I ) =PZ( I ) 

RX ( I ) =RX( I )+OTOD*( HDTOM*FX ( I )+VX ( I ) ) 

RY ( I ) =RY( I )+DTOD*( HDTOM*FY( I ) + VY ( I ) ) 

RZ ( I ) = RZ ( I ) +DTOQ*( HDTOM*FZ ( I ) +VZ ( I ) ) 

245 CONTINUE 
CALL STEP 
EM AX = 0 . 0 
FM AX=0 . 0 
T I ME=T IME+DT 
NT = NT + 1 

IF(LCUTd) .GT.O) GO TO 265 
1=1 

VSS = VX ( I ) 

VX( I )=VSS+HDTOMB*FX( I ) 

RX ( I ) = RXK ( I ) +( VX ( I ) + VSS ) *HDTOD 
VSS=VY ( I ) 

VY ( I ) =VSS+HDT0M8*FY( I ) 

RY ( I ) =RYK ( I )+(VY( I ) +VSS ) *HDTCD 
VSS=VZ( I ) 

VZ( I )=VSS+HDTOMB*FZ< I ) 

RZ ( I ) = RZK ( I ) + ( VZ ( I )+VSS )*HDTOO 

PKE ( I ) =VX( I )* V X ( I ) * VY ( I ) *VY ( I ) +VZ( I ) *VZ( I ) 

EMAX1=PKE ( I ) 

FMAX 1 = FX { I )*FX( I ) + FY ( I )*FYU)+FZ( I )*FZ( I ) 
DTE1=TEFAC*SQRT{ 1.0/FMAX1 ) 

DTF1=SQRT{ TFACB/FMAX1 ) 

260 FX ( I ) = 0 . 0 
FY( I ) = 0. 0 
FZ ( I )=0.0 
265 DO 280 1=2, LD 

IF(LCUT( I ) .GT. 0 ) GO TO 280 
V SS = VX( I ) 

VX( I ) = VSS+HDTOM* FX ( I ) 

RX Cl) =RXK ( I ) + ( VX ( I ) + VSS)*HDTCD 
VSS = VY ( I ) 

VY ( I )=VSS+HDTOM*FY ( I ) 

RY( I ) =RYK ( I ) + ( VY( I ) + VS S ) *HDT OD 
VSS=VZ( I ) 

VZ (I ) = V S S + HDTO M * p Z ( I ) 

RZ( I )=RZK( I ) + { VZ l I ) +VSS)*HDTOD 
PKE ( I ) = V X C I)*VX( I ) +VY ( I )*VY( I ) + VZ ( I ) *V Z ( I ) 
275 F2 = FX ( I ) *FX ( I )+FY{ I ) * FY ( I ) +FZ ( I ) *F Z ( I ) 

FX { I ) = 0. 0 
FY ( I ) = 0 . 0 
FZ ( I )=0.0 

I F ( F2.GT.FMAX) FMA X=F2 
IF (PKE ( I ) .GT .EMAX) EMAX=PKE(II 
280 CONTINUE 
DTL=DT 

CT I ME=0. 01*1 TIME (XX) -START 
DTE=TE FAC* SORT (1 .O/EMAX) 

DTF=SQRT (TFAC/FMAX) 

I F ( EMAX1 . GT . EM AX ) EMAX=EMA X 1 
DT=DTE 1 

IF(DT.GT.DTFl) DT = DTF1 
IF (DT.GT .DTE ) DT=DTE 
IF(DT.GT.DTF) DT = DTF 
IFIDT.GT.1.0E-14) DT=1.0E-14 
300 I F ( ISHUT.EQ.-l ) GO TO 400 
310 IF(NS-NT) 400,400,320 
320 DO 325 1=1, LL 

VX ( I ) = V F AC*VX ( I ) 

VY ( I ) =VFAC*VY( I ) 
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325 VZ ( I )=VFAC*VZ( I ) 

GO TO 800 
C 

THE PRINT SUBROUTINE PLACES A HEADING OF PERTINENT INFORMA- 
TION AT THE TOP OF EACH TIMESTEP PRINTOUT. 

POTENTIAL ENERGY AND LOCAL POTENTIAL ENERGY FOR EACH ATOM 
ARE CALCULATED BASED ON THE NEW POSITIONS. SUMMATIONS OF 
TOTAL POTENTIAL AND KINETIC ENERGY FOR THE LATTICE ARE PER- 
FORMED. DX, DY, AND DZ KEEP TRACK OF MOTION RELATIVE TO THE 
INITIAL POSITION AT TIME ZERO FOR EACH ATOM. 

400 CALL PRINT 
C 

410 T POT = 0.0 

DO 450 1=1, LL 
PP E ( I ) =0 .0 
450 PT E ( I ) =0 .0 
CALL ENERGY 
PKE( 1 )=HBMAS*PKE ( 1 } 

TPKE=PKE( 1 ) 

PT E ( 1) =PKE ( 1)+PPE( 1) 

DO 620 1=2, LL 

PKE ( I ) =HTMAS*PKE ( I ) 

TPKE=TPKE+PKE( I ) 

620 PTE( I ) =PKE( I ) +PP E ( I) 

TE=TP0T+8 I ND 
WRITE ( 6,9660) 

DTEST=(RY( lJ-RYI ( 1) )**2 
IF (DTEST.CT. 0.01) DT EST= 0.01 
I F ( T POT . LE.TPOTL) GO TO 700 
ERAT =T PKE-RLL 
700 DO 750 1=1, LD 

DX ( I )=RX( I ) -RX I ( I ) 

DY ( I ) = RY ( I ) -RY I ( I ) 

DZ( I )=RZ( I J-RZI ( I ) 

IF (DX(I)**2.GE.DTEST) GO TO 720 

IF <DY(I)**2.GE.DTEST) GO TO 720 

IF (DZ( I)**2.GE.DTEST) GO TO 720 

GO TO 750 
C 

THIS SECTION PRINTS THE RELATIVE MOTION, VELOCITY, AND 
ENERGY OF EACH ATOM, FOR EVERY TIMESTEP SO DESIGNATED: IE, 

EVERY ND'TH TIMESTEP, BEGINNING WITH U NS AND ENDING WITH 
#NTT . 

720 WRITE ( 6,9670) I , DX ( I ) , D Y ( I ) , D Z ( I ) , V X ( I ) , VY ( I ) , 

1VZ( I ) , PKE ( I ) ,PPE ( I ) , PT E ( I) 

750 CONTINUE 

WRITE ( 6 ,9640) T PK E , T POT , T E, RO EL 
WRITE ( 6,9650) NPAGE 
NPAGE=NPAGE+1 
TPOTL=TPOT 

IF ( NT-NTT) 760,950,950 
760 DO 780 1=1, LL 

VX(I ) =VFAC*VX ( I ) 

VY ( I )=VFAC*VY( I ) 

VZ( I )=VFAC*VZ( I ) 

780 CONTINUE 

I F ( ISHUT.EQ.— 1) GO TO 950 
790 NS=NS+ND 

800 IF (NDtC.EQ.10) GO 1C 810 
GO TO 100 
810 DDT I =0 • 1 *DDT I 
DT I =DT I +DDT I 
NDEC=0 
GO TO 100 
950 CONTINUE 
C 

THIS SECTION PRINTS OUT X, Y, AND Z COORDINATES AND BINDING 
ENERGIES OF EACH ATOM IN THE CRYSTAL AT THE END OF THE 
PROGRAM. ALSO, DATA CAPOS ARE PRINTED WITH X,Y,Z COORDINATE 
OF EACH ATOM IN THE CRYSTAL FOR USE IN THE DYNAMIC PROGRAM. 
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955 WRITE ( 6,9620) IH2,NT 

WRITE (7,9690) LL , D1X , D 1Y , D 1Z , NVAC 
DO 965 1=1, LL, 3 
K= I + 1 
J= I +2 

WRITE (7,9691) RX( I ) , RY ( I ) , RZ ( I ) , RX ( K) , RY ( K ) , RZ ( K ) , RX 
1 RZ ( J ) 

965 WRITE ( 6 ,9630) I , RX ( I ) , RY ( I ) , RZ ( I ) , P PE ( I ) , K , R X ( K ) , 
1RY(K) ,RZ(K) ,PPE(K) ,J,RX(J) , RY ( J ) ,RZ(J) , PP E ( J ) 

WRITE ( 6,9640) TPKE , TPOT , TE , ROEL 
WRITE ( 6,9650) NPAGE 
1000 IF C I SHUT ) 9999,5,5 
9999 STOP 
END 



SUBROUTINE CROSYM 
C 

SOLVES M SIMULTANEOUS EQUATIONS BY THE METHOD OF CROUT 
THIS SUBROUTINE FITS THE BEST CUBIC BETWEEN THE REPULSIVE 
AND ATTRACTIVE PARTS OF THE POTENTIAL. 

C 

COMMON/COMA/ A(4,5),MCR0 
M=MCRO 
N=M+ 1 
11 = 1 

100 13=11 

SUM= AB S ( A ( 11,11) ) 

DO 120 1=11, M 

I F ( SUM- ABS ( A ( I ,11))) 110,120,120 
110 13=1 

SUM= ABS ( A ( I, ID ) 

120 CONTINUE 

I F ( 1 3- I 1 ) 130, 150, 130 

130 DO 140 J = 1 , N 
SUM=— A ( I 1 , J) 

A( II, J )=A( 13, J ) 

140 A ( 13 , J )=SUM 
150 13=11+1 

DO 160 1=13, M 

160 A( I , 11 )=A( I , II )/A( II, I 1) 

170 J 2= I 1-1 
13=11+1 

I F ( J 2 ) 180,200,180 

180 DO 190 J=I 3,N 
DO 190 1=1, J2 

190 A( H , J )=A( II , J ) — A( II , I )*A( I , J) 

I F ( 1 1-M ) 200,220,200 
200 J2 = 1 1 
11=1 1+1 
DO 210 1=1 1 , M 
DO 210 J = 1 , J 2 

210 A( I , 1 1 ) = A ( I, Il)-A( I , J ) * A ( J , II) 

IF(Il-M) 100,170,100 
220 DO 240 1=1, M 
J2=M-I 
1 3 = J 2+ 1 

A ( I 3 , N ) = A ( 1 3 , N ) / A ( 13, 13) 

I F ( J 2) 230,250,230 

230 DO 240 J=1,J2 

240 A( J,N)=A( J,N)-A( I3,N)*A{ J, 13) 

250 RETURN 
END 
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SUBROUTINE B100 

THIS IS A LATTICE GENERATOR THE THE BCC (100) ORIENTATION. 
THE CRYSTAL IS DEVELOPED IN THE ORDER, X FOLLOWED BY Z, 
FOLLOWED BY Y. 

IT CONTAINS A NONSTANDARD USE OF THE SURFACE RELAXATION 
PARAMETER. 

C 

C0MM0N/CCM1/RX (1000) ,RY( 1000 ), RZ( 1000) ,LCUT( 1000) , 
ILL , LD , I TYP E , NVAC 

COMMON /C0M4/ IX, IY, I Z , SCX , SC Y , SC Z , I DE EP , D1 X , D 1 Y , D1 Z , 

1 I V AC X , I VAC Y , IVACZ 
DIMENSION YL AX ( 20 ) 

DATA Y LAX/ 20*0 . 0/ 

YLAX ( 1 ) =-G .20 
YLAX ( 2 ) =-0.03 
SCX=1.0 
SC Y = 1 . 0 
SCZ=1. 0 
M = 2 
JT = 0 
Y=-S CY 

DO 60 J=l, IY 
Y=Y+SC Y 
KT = 0 
Z=-SCZ 

DO 59 K=1 , I Z 
Z=Z+SCZ 
IT = 0 
X=-SCX 

DO 58 1 = 1, IX 

IF( IT-( IT/2)*2) 21,11,21 

11 IF(JT-(JT/2)*2) 57,12,57 

12 IF (KT- (KT/2) *2 ) 57,30,57 

21 I F ( J T- ( J T / 2 ) * 2 ) 22,57,22 

22 IF (KT- ( KT/2 ) *2 ) 30,57,30 
30 RX(M)=X 

RY ( M ) =Y+YLAX ( J ) 

RZ ( M ) = Z 
M=M+ 1 

IF ( IT .NE. IVACX) GO TO 57 
IF ( JT .NE. IVACY ) GO TO 57 
IF (KT.NE. IVACZ) GO TO 57 
NVAC=M- 1 

57 I T= I T + 1 

58 CONTINUE 
KT = KT+ 1 

59 CONTINUE 

JT = JT + 1 

IF(IDEEP-JT) 60,110,60 

60 CONTINUE 
LL=M-1 

100 RETURN 
110 LD=M- 1 

GO TO 60 
END 



SUBROUTINE PLACE 
C 

THIS SUBROUTINE LOCATES A VACANCY, INTERSTITIAL, OR REPLACE- 
MENT IMPURITY IN THE LATTICE. 

C 

COMMON/COM 1/RX ( 1000) ,RY( 1000 ), RZ ( 1000) ,LCUT( 1000) , 

1 LL , LD , I TYP E , NVAC 

COMMON /C0M4/ IX , I Y, IZ , SCX , SCY ,-S C Z , I DE EP , D1 X , D1Y , D1 Z , 

1 IVACX , IVACY, IVACZ 
GO TO (10,20,30,40), ITYPE 
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10 LCUT(NVAC) = 1 
LCUT(l) = 1 
RX(1) =0. 0 
RY { 1 ) = 0 . 0 
RZ ( 1 ) = 0 . 0 
GO TO 50 

20 RX ( 1 )=RX(NVAC ) +D1X 
RY(1 )=RY (NVAC) +D1Y 
RZ ( 1 ) = RZ (NVAC) + D1Z 
GO TO 50 

30 LCUT(NVAC) = 1 
R X ( 1 ) = RX (NVAC) 

RY ( 1 ) = P Y (NVAC ) 

RZ ( 1 ) = RZ(NVAC) 

GO TO 50 

40 RX (1 )=RX (NVAC) +D1X 
RY(1)=RY(NVAC)+D1Y 
RZ ( 1 ) = RZ (NVAC) +D1Z 

50 RETURN 
END 



SUBROUTINE STEP 
C 

THIS SUBROUTINE DOES THE DYNAMICS FOR ONE TIMESTEP. 

THE FIRST HALF DOES THE DYNAMICS FOR ATOM #1; THE SECOND 
HALF FOR ALL OTHERS. 

C 

COMM ON/ C OM 1 /RX ( 100 0) ,RY( 1 000 ) , R Z ( 1 000 ) , LCUT( 1000) , 
1LL,LD, ITYPE,NVAC 

COMMON /C0M5/R0E , R0E2, ROEM, EXA,EXB, PEXA,PEXB» FXA,PFXA, 
HQtTSAVEt BSAVE 

COMMON /C0M6/FX ( 1000) , FY ( 1 000 ) , FZ ( 1 000 ) , PAC , P FPTC , FM 
COMMON /C OM 8/ROE A , ROE B , ROEC , R0EC2 ,C PO ,C PI ,CP2,CP3, 
1CF0, CF1 ,CF2,CGD1 , CG02, CGB1 , CGB2 , CGF 1 , C GF 2 
IF (IQ-1) 100,101,102 

100 I P = 2 

GO TO 200 

101 IP = 1 

GO TO 200 

102 1=1 

I P = 2 

105 DO 195 J = I P , L L 

I F ( LCUT ( J ) ) 195,110,195 
110 DR X=RX ( J ) -RX ( I ) 

IF(DRX) 113,117,117 
113 IF ( DRX+ROE ) 195,195,120 
117 I F LDRX-ROE ) 120,195,195 
120 DRY=RY ( J ) -RY ( I ) 

I F ( DRY ) 123,127,127 
123 IF ( DRY+ROE ) 195,195,130 
127 I F ( DRY-ROE ) 130,195, 195 

130 DRZ=RZ( J)-RZ(I) 

IF(DRZ) 133,137,137 
133 IF { DRZ+ROE ) 195,1 Q 5,140 
137 I F ( DRZ-ROE ) 140,195,195 
140 DIST=DRX*DRX+DRY*DRY+DRZ*DRZ 
I F ( D I S T-RO E2 ) 150,195, 195 
150 DI ST = SQR1 (DJ ST) 

160 IF ( DI ST- ROEM ) 162,162,165 

162 F0RCE=2XP ( PFXA+P EX 8*DI ST ) 

GO TO 180 
165 DFF= RO E-D I ST 

I F ( DFF-1 .0E-10 ) 195, 195, 167 

167 FORC E= ( E XP ( D AC + PEXB*DI ST ) -PFPTC ) /OFF 
180 IF ( FM- FORCE ) 190,190,195 

190 FOD=FORCE/DIST 
FA=FOD*DRX 
FX ( J ) = FX ( J ) +F A 
FX ( I ) =FX ( I ) -F A 
FA=FOD*DRY 
FY ( J ) =FY ( J ) +FA 
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FY ( I ) =FY ( I ) -FA 
FA=FOD*DRZ 
FZ( J ) = FZ( J J+FA 
FZ t I ) = F Z ( I } — FA 
195 CONTINUE 

200 DO 300 I = I P , LD 

I F ( LCUT ( I ) ) 300,205,300 
205 IP=I+1 

DO 295 J=I P , LL 
I F ( LCUT ( J ) ) 295,210,295 
210 DRX=RX ( J )-RX ( I ) 

IF(DRX) 213,217,217 
213 I F ( DRX+ROFC ) 295,295,220 
217 IF ( DRX-ROEC) 220,295,295 
220 DRY=RY ( J )-RY { I ) 

I F ( DRY ) 223,227,227 
223 I F ( DRY+ROEC ) 295,295,230 
227 I F ( DRY-ROEC) 230,295,295 
230 DRZ = RZ ( J )-RZ ( I ) 

IF(DRZ) 233,237,237 
233 IF (DRZ+RGEC) 295,295,240 
237 IF ( DRZ-ROE C ) 240,295,295 
240 DIST=DRX*DRX+DPY*DRY+DRZ*DRZ 
I F ( D I S T-R0EC2 ) 250,295,295 
250 DIST=SQRT(DI ST) 

I F { D I ST-Ru EA ) 260,255, 255 
255 IF ( D I ST-ROEB) 265,270,270 
260 FORCE=EXP(FXA+EXB*DIST) 

GO TO 280 

265 FORCE=DIST*(DI ST*CF2+CF1 ) +CFO 
GO TO 200 

270 F0RCE=EXP(CGF1+CGB1*DIST )-EXP(CGF2+CGB2*DI ST ) 
280 IFtABS(FORCE) . LE.FM) GO TO 295 
FOD = FORC E/DI ST 
FA=FOD*DRX 
FX(J)=FX(J)+FA 
FX ( I ) = FX ( I ) -FA 
FA = F OD^-'DRY 
FY( J ) =F Y( J ) +FA 
FY ( I ) = FY ( I ) -F A 
FA=FOD*DRZ 
FZ( J)=FZ( JH-FA 
FZ ( I ) = FZ ( I )-F A 
295 CONTINUE 
300 CONTINUE 
RETURN 
END 



SUBROUTINE ENERGY 
C 

THIS SUBROUTINE CALCULATES THE MUTUAL POTENTIAL ENERGIES. 
THE FIRST HALF DOES THE DYNAMICS FOR ATOM #1; THE SECOND 
HALF FOR ALL OTHERS. 

C 

COMMOM/COMl/RXI 1000) ,RY(1000) ,RZ(1000) , LCUT ( 1000) , 

1 LL , LD , ITYPE.NVAC 

COMM GN /CQM5/ ROE, RO E2 ,ROEM, EXA, EXB, PEXA, PEXB , FXA, PFXA, 
1IQ,TSAVE,BSAVE 

COMMON /C0M7/P P TC, T POT, P»E( 1000 ) , TL P E , R OEL , R0 C L2,NEW 
COMMON /COM 8/ROE A, ROEB,ROEC , R0EC2 ,C PO ,CP1 , C P2 , CP3 , 
1CF0, CF 1,CF2,CGD1 »CGD2»CG31 ,CG32,CGF 1 ,CGF2 
IF (IQ-1) 100,101,102 

100 I P=2 

GO TO 200 

101 IP=1 

GO TO 200 

102 1=1 

I P=2 

105 DO 595 J = 1 P , LL 

I F ( LCUT ( J ) ) 595, 510, 595 
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510 DRX= RX ( J ) -RX ( I ) 

IF(DRX) 513,517,517 
513 I F ( DRX + ROE ) 595,595, 520 
517 I F ( DRX-RQE ) 520,595,595 
520 DR Y=RY ( J ) — RY ( I ) 

I F ( DRY ) 523,527,527 
523 I F ( DRY +ROE ) 595,595,530 
527 I F ( DRY-ROt ) 530,595,595 
530 DRZ^RZ ( J ) — RZ ( I ) 

IF(DRZ) 533,537,537 
533 I F ( DRZ+ROE ) 595,595,540 
537 I F ( DRZ-RGE ) 540,595,595 
540 D 1 ST=DRX*DRX+DRY*DRY+DRZ*DRZ 
I F ( D I ST-RO E2 ) 550, 595,595 

550 DIST=SQRT (DIST) 

560 POT=EXP(PEXA+PEXB*DIST)-PPTC 
580 TPOT=TPOT+POT 

PPE ( 1 ) =PPE ( I ) + BS AV E*PO T 
PPF ( J ) =PPE ( J )+TSAVE*POT 
595 CONTINUE 
600 CONTINUE 

200 DO 300 I=IP,LD 

IF (LCUT( I ) ) 300,205, 300 
205 IP=I+1 

DO 295 J = I P , L L 
IF (LCUT ( J ) ) 295,210,295 
210 DRX=RX ( J ) -RX ( I ) 

IF(DRX) 213,217,217 
213 IF ( DRX +ROEC ) 295,295,220 
217 I F ( DRX-RGEC) 220,295,295 
220 DRY= RY ( J ) — R Y ( I ) 

I F ( DRY ) 223,227,227 
223 I F ( DR Y+ROEC ) 295,295,230 
227 IF ( DRY-ROEC) 230,295,295 
230 DRZ=RZ ( J )— RZ ( I J 

IF(DRZ) 233,237,237 
233 I F ( DRZ +ROEC J 295,295,240 
237 IF(DRZ-ROEC) 240,295,295 
240 DI ST = DRX- : DRX+DRY*DRY+DRZ*D RZ 
IF ( DIST-R0EC2) 250,255,295 
250 DI ST=SQRT (DI ST ) 

IF(DIST-ROEA) 260,255,255 
255 IF ( DIST-ROEB) 265,270,270 
260 POT=EXP ( EXA+EXB*DI ST ) 

GO TO 280 

265 POT=DIST*( DIST*( DI ST*C P3 +C P 2 ) +C P 1 ) +CPO 
GO -TO 280 

270 P0T=EXP(CGD1+CGB1*DIST)-EXP(CGD2+CGB2*DI ST) 
280 TPOT=TPOT+POT 
SAVE = 0 . 5* POT 
PP E ( I ) =PPE ( I ) +SA VE 
PPE(J)=PPE(J )+SAVE 
295 CONTINUE 
300 CONTINUE 
RETURN 
END 



SUBROUTINE PRINT 
C 

THIS SUBROUTINE PRINTS THE HEADING OF ALL PERTINENT INFORMA 
TION AT THE TOP OF EACH TIMESTEP PRINTOUT. 

C 

COMMON/C CM !/RX(1000) ,RY( 100 0),RZ(1000) ,LCUT( 1000) , 

ILL , L D, ITYPE,NVAC 

C0MM0N/C0M2/ I HI ( 20 ) , I H2 ( 8 ) , I HS < 10 ) , I HB ( 6 ) , IHT(6) , 

1 TARGET (4) , TMAS , BULLET ( 4 ) , B M AS-, P LAN E , T E MP , T HE RM 
COMMON /C0M3/RX I ( 1000) ,RYI( 1000) ,RZI ( 1000) ,CVR,EVR, 

1NT , T I ME , OT , DT I , I LAY 

COMMON /COM 4/ IX,IY,IZ,SCX,SCY,SCZ,IDEEP,D1X,D1Y,D1Z, 
1IVACX, 1VACY, I VAC Z 
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PRIMARY ENERGY 
, I 2, 3H X , 12, 3H 
14/) 

PRIMARY ENERGY 
, 12, 3H X , 12, 3H 



X , 12, 3H 



X\ 12,3 
FROM 



, 12 ,3H X , 12, 3H 
=,F5.2,5H, Y =, 



COMM ON/C OM5/ROE , R0F2,R0EM, EX A, EXB, P EXA , P EXB , F XA , P F XA , 
1IQ,TSAVE ,BSAVE 

COMMON/ COM8/ROFA,ROEB,ROEC, ROEC 2 ,CP0»CPL,CP2,CP3» 
1CF0,CF1,CF2,CGDI ,CGD2, CGB1 , CGB2 , CGF1, CGF2 

9710 FORMAT (40X, 10A4, /, 28X, 20A4, /) 

9720 FORMAT (9H TARGET -,4A4, 10HPRIMARY - , 4A4 , IX , 14HLATTI CE 
1 UNIT = , F 7 . 4 , 4H ANG) 

9730 FORMAT (4X,6HMASS = , F7 . 2 , 13 X , 6HMA SS = , F7 . 2 , 9X , 14H L ATT I C 
IE TEMP = F5 .2 ,7H DEG K,,18H THERMAL CUTOFF =,F5.2,3H E 
IV/) 

9740 FORMAT (2H (,A4,8H) PLANE,, 18H 
1 F5 . 2 , 21HK EV , CRYSTAL SIZE ( 

1 ),, 4 X , 16H VACANCY IN SITE , 

9741 FORMAT ( 2H (,A4,8H) PLANE, , 18H 
1 F5.3,21HKEV, CRYSYAL SIZE ( 

1H ),, 4X, • INTERSTI TI AL ( • , 2 ( F5. 2 , ' , • ) , F5.2 , • ) 

1SITE • ,14/ ) 

9742 FORMAT ( 2H (,A4,3H) PLANE, ,18H PRIMARY ENERGY =, 

1 F 5 . 2 , 2 1HKE V , CRYSTAL SIZE ( ,I2,3H X 

1 ),, 4X , 20HREPLACEMENT IN SITE , 14/) 

9750 FORMAT ( 30H PRIMARY START POINT ( LU ) X 

1F5.2,5H, Z =,F5. 2, 5X, 13, • LAYERS ARE FREE TO MOVE', 

1 10X , 4H I Q =,12/ ) 

9760 FORMAT ( 12H POTENTIAL , 6 A4 , 3X , 5H PE X A= , F9 . 5 , 2X , 5HP EXB = , 
1F9.5,2X,5HPFXA=,F9.5) 

9765 FORMAT (12X,6A4,3X, 5HEXA =, F9 . 5 , 2X, 5HEXB = , F9 . 5 , 2X , 5HFX 
1A = , F9 . 5/ ) 

9770 FORMAT! • WHEN • , F8. 4, ' < R <*,F8.4,' 

1TIAL PARAMETERS ARE',//,' CPO =',F10.3, 

IF 1 0. 3 , ' , CP2 = ' , F10.3,', CP3 =',F10.3,/,' 

1E10.3,', CF1 = ' , ElO. 3, ' , CF2 =',E10.3,//) 

9780 FORMAT! • CUT-OFF AT',F5.2,', WHEN R > ',F6.3,' L'J, MOR 
1SE POTENTIAL PARAMETERS ARE', 8A4,//,10X,' CGD1 =' , 
1F8.4,', CGD2 = ' , F8.4, ' , CGB1 =',F8.4,', CGB2 =',F8.4, 
1', CGF 1 = ' , F 8 . 4 , ' , CGF2 =',F8.4,//) 

9790 FORMAT! 10H TIMESTEP , 1 4 , 22X , 6HDTI = , F5.4, 5H LU, 

1 , 22H ELAPSED TIME (SEC) =, ElO. 4,', NEXT TIMESTEP IS 
1=' , E 1 0 . 4/ ) 

WRITE ( 0,9710) I HS , I HI 
WRITE ( 6,9720) T ARGET , BU L L ET , CVR 
WRITE ( 6,9730) TMA S , BMAS , TEMP , TH ERM 
GO TO (401,402,403,402), I TYPE 

401 WRITE ( 6,9740) P LAN E , E VR , I X , I Y , I Z , NV AC 
GO TO 405 

402 WRITE ( 6,9741) 

GO TO 405 

403 WRITE ( 6,9742) P LAME , E VR , I X , I Y , I Z , NVAC 

405 WRITE ( 6 ,9750) R X I ( 1 ) , RY I ( 1 ) , R Z I ( 1 ) , I L A Y , I Q 
( 6,9760) IHB,PEXA,PEXB,PFXA 

( 6,9765) IHT ,EXA,EX3,FXA 

( 6,9770) ROEA, ROEB , CPO , CPI , CP2, CP3 , CFO,CF 1 , CF2 
( 6,9780) ROEC , ROEB , IH2 ,CGD1 ,CGD2 ,CGB1 ,CGE2, 



THE MATCHING POTEN 
CPI =• 
CFO =' 



PLANE , EVR , IX, IY,IZ,D1X,D1Y,D1Z,NVAC 



WRITE 
WRITE 
WRITE 
WRITE 
WRITE 

1CGF1 , CGF2 

WRITE ( 6,9790) 
RETURN 
END 



NT , DT I , T I ME, DT 



BLOCK DATA 

DIMENSIONING OF VARIABLES USED IN COMMON 

COMMON/ C0M1/RX ( 1000 ), RY! 1000 ), R Z ( 1 000 ) ,LCUT( 1000) , 

ILL, LD, ITYPE, NVAC 

DATA RX/1000*0.0/»RY/1000*0.0/,RZ/1000*0.0/,LCUT/1000* 
COMMON/ C0M3/RX I (1000), RY I ( 1000 ) , RZ I ( 1000) , CVR, EVR, 

1NT, TIME, DT, DTI , I LA Y 

DATA RXI/ 100 0*0.0/ ,RY! / 1000*0. 0/,RZI/1000*0.0/ 
C0MM0N/C0M6/FX ( 1000 ), FY ( 1000 ), FZ ( 1 000 ) , PAC , P F°TC , FM 
DATA F X/ 1000*0. 0/,FY/1000* 0.0/ , FZ/ 1000*0.0/ 

END 

//G0.FT06F001 DD SPACE=(CYL, (10,2) ,RLSE) 

//GO . SYS I N DD * 

CRYSTAL-1968 MODIFIED TO DEAL WITH VACANCIES AND I NT5RS.TI TI 
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( GIRIFALCO — WEIZER POTENTIAL ) .99060 

ARGON 39.948 9.33 -5.60 

TUNGSTEN 183.86 11.30 -7.50 

BODY CENTERED CUBIC, (100) ORIENTATION 
100 5 4 51 



1.41160 3.03200 3.4 
ARGON-TUNGSTEN 
TUNGSTEN-TUNGSTEN 
100 10 10 10 
0.7 -0.7 0.0 5 1 
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